/*************************************************************
** THROBIT1: Threshold observability bivariate probit model
** Monte Carlo Comparison of Throbit with "Boolean" Probit
** (c) 2003 by Sanford C. Gordon
**
**
** Department of Politics
** new York University
** Purpose: Determine sample properties of partially observed probit with thresholds.
** output file: mc1.out
*************************************************************/

rndseed 33333;

closeall; 
cls; 

/* Call to the Maxlik library */
library maxlik;
/* Resets global variables to default values */
maxset;

/* Global Parameter Settings */

_max_algorithm = 2;
	                 	@ Optimization Method           @
                       	@ 1 = Steepest Descent          @
                       	@ 2 = BFGS                      @
                       	@ 3 = DFP                       @
                       	@ 4 = Newton-Raphson            @
                       	@ 5 = BHHH                      @
                       	@ 6 = PRCG                      @

_max_CovPar = 1;
						@ Covariance Matrix Method		@
						@ 0 = Not Computed			    @
						@ 1 = Numerical Hessian			@
						@ 2 = Quasi-maximum likelihood  @

_max_LineSearch = 2;
						@ Line Search Method			@
						@ 1 = unit step length			@
						@ 2 = cubic or quadratic		@
						@ 3 = step halving			    @
						@ 4 = Brent's method			@
						@ 5 = BHHH method				@
__output=0;
/*Create symbolic infinity*/
infty=2^10000;

proc throbiti(parms,dta);
    local obs,y,y1,y2,x1,x2,beta1,beta2,tau1,tau2,idx1,idx2,p0,py1,py2,pya,llike;
    /* Data */
    obs = rows(dta);
    y   = dta[.,1];
    y1  = dta[.,2];
    y2  = dta[.,3];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    /* parameters */
    beta1 = parms[1:3];
    beta2 = parms[4:6];
    tau1  = exp(parms[7]);
    tau2  = exp(parms[8]);
    idx1= x1*beta1;
    idx2=x2*beta2;
    p0 = cdfn(-idx1).*cdfn(-idx2);
    py1 = cdfn(idx1-tau1).*cdfn(tau2-idx2);
    py2 = cdfn(idx2-tau2).*cdfn(tau1-idx1);
    pya = 1-p0-py1-py2;
    llike = (1-y).*ln(p0)+y.*y1.*ln(py1)+y.*y2.*ln(py2)+y.*(1-y1).*(1-y2).*ln(pya);
    retp(llike);
endp;

proc boolprob(parms,dta);
    local y,llike,obs,x1,x2,beta1,beta2,p1,p2,yesprob,noprob;
    obs=rows(dta);
    y=dta[.,1];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    beta1=parms[1:3];
    beta2=parms[4:6];
    p1 = cdfn(x1*beta1);
    p2 = cdfn(x2*beta2);
    yesprob = 1-(1-p1).*(1-p2);
    noprob  = 1-yesprob;
    llike = y.*ln(yesprob)+(1-y).*ln(noprob);
    retp(llike);
endp;

/* Parameter Values Fixed for all simulations */

beta1=1|1|1;
beta2=1|1|1;
tau1 = 2;
tau2 = 2;
rho=0;
mciter=1000;

/**************************************************************************************
** First Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 100  **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=100;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls;
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));
    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = (thetath-truthth)*(thetath-truthth)';
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp100 = meanc(packr(serrbp));
seth100 = meanc(packr(serrth));

rmsebp = sqrt(diag(msematbp))/mciter;
rmseth = sqrt(diag(msematth))/mciter;

/* Print table */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meanth = meanc(resultth);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
format 200,2;
output file = "c:/gauss50/research/throbit/mc1.out" reset;
print"\\begin{table}";
print"\\caption{\\label{mc1} Monte Carlo Comparison of Boolean Probit and Throbit}";
print"\\begin{center}";
print"\\begin{tabular}{lcccccccc}";
print"\\hline \\hline";
print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{13}$&$\\beta_{20}$&$\\beta_{2c}$&$\\beta_{23}$&$\\ln(\\tau_1)$&$\\ln(\\tau_2)$ \\\\ \\hline";

print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=100}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Proportion successful convergence, Throbit: " ""$+ftocv(rows(packr(serrth))/rows(serrth),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diag.out" reset;
print "I have finished with N=100, moving on to N=500";
output off;



/**************************************************************************************
** Second Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 500 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=500;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));
    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = (thetath-truthth)*(thetath-truthth)';
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp500 = meanc(packr(serrbp));
seth500 = meanc(packr(serrth));

rmsebp = sqrt(diag(msematbp))/mciter;
rmseth = sqrt(diag(msematth))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meanth = meanc(resultth);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=500}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Proportion successful convergence, Throbit: " ""$+ftocv(rows(packr(serrth))/rows(serrth),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diag.out" on;
print "I have finished with N=500, moving on to N=1000";
output off;



/**************************************************************************************
** Third Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 1000 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=1000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));
    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = (thetath-truthth)*(thetath-truthth)';
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp1000 = meanc(packr(serrbp));
seth1000 = meanc(packr(serrth));

rmsebp = sqrt(diag(msematbp))/mciter;
rmseth = sqrt(diag(msematth))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meanth = meanc(resultth);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Proportion successful convergence, Throbit: " ""$+ftocv(rows(packr(serrth))/rows(serrth),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diag.out" on;
print "I have finished with N=1000, moving on to N=5000";
output off;

/**************************************************************************************
** Fourth Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 5000 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=5000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));
    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = (thetath-truthth)*(thetath-truthth)';
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp5000 = meanc(packr(serrbp));
seth5000 = meanc(packr(serrth));

rmsebp = sqrt(diag(msematbp))/mciter;
rmseth = sqrt(diag(msematth))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meanth = meanc(resultth);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";


N = 1000000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1  = ones(N,1)~xcommon~x12;
x2  = ones(N,1)~xcommon~x22;
epsilon1 = rndn(N,1);
epsilon2 = rndn(N,1);
zstar1   = x1*beta1+epsilon1;
zstar2   = x2*beta2+epsilon2;
y  = maxc((zstar1~zstar2)').>0;
y1 = (zstar1.>tau1).*(zstar2.<tau2);
y2 = (zstar2.>tau2).*(zstar1.<tau1);
yprop  = meanc(y-y1-y2);
y1prop = meanc(y1);
y2prop = meanc(y2);

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=5,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Proportion successful convergence, Throbit: " ""$+ftocv(rows(packr(serrth))/rows(serrth),0,2)$+"}\\\\";
print " \\hline \\hline ";
print "\\multicolumn{9}{l}{Approximate proportion ambiguous success: " ""$+ftocv(yprop,0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Approximate proportion success for first reason: " ""$+ftocv(y1prop,0,2)$+"}\\\\";
print "\\multicolumn{9}{l}{Approximate proportion success for second reason: " ""$+ftocv(y2prop,0,2)$+"}\\\\";
print "\\end{tabular}";
print "\\end{center}";
print "\\end{table}";
output off;

output file = "c:/gauss50/research/throbit/diag.out" on;
print "I have finished with N=5000.";
print "Mean standard errors, N=100";
print "Boolean probit:" sebp100';
print "throbit:" seth100';

print "Mean standard errors, N=500";
print "Boolean probit:" sebp500';
print "throbit:" seth500';

print "Mean standard errors, N=1000";
print "Boolean probit:" sebp1000';
print "throbit:" seth1000';

print "Mean standard errors, N=5000";
print "Boolean probit:" sebp5000';
print "throbit:" seth5000';
output off;





